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Abstract. We study the non equilibrium time evolution of an integrable field theory 
in 1 + 1 dimensions after a sudden variation of a global parameter of the Hamiltonian. 
For a class of quenches defined in the text, we compute the long times limit of the one 
point function of a local operator as a series of form factors. Even if some subtleties 
force us to handle this result with care, there is a strong evidence that for long times 
the expectation value of any local operator can be described by a generalized Gibbs 
ensemble with a different effective temperature for each eigenmode. 
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1. Introduction 

In recent years, the unitary time evolution of an extended quantum system has attracted 
a lot of attention. The interest of the community on this topic was spurred by novel 
experiments with cold atoms that have shown the coherent evolution of a quantum 
system in a laboratory [IHI]. One of the simplest and most studied setting is the so 
called (sudden) quantum quench. In this situation, one initially considers an extended 
quantum system prepared in a pure state: this can be regarded as the zero temperature 
ground state of an Hamiltonian At t=0 the Hamiltonian is suddenly changed from 
"Ho to % and then the system is let to evolve unitarily according to the new Hamiltonian 
"H, without any coupling to the environment. Sure enough, there is a transient regime for 
such an abrupt change of the Hamiltonian but the main interest is in what happens in the 
long time limit. A relevant issue is whether the expectation values of local operators at 
t = oo can be derived by standard thermodynamical ensembles (e. g. Gibbs ensemble): 
if this is the case, it is said that the system thermalizes. The usual understanding is 
that thermalization should generally occur in extended systems but this notion has been 
recently challenged by experiments done on one-dimensional quasi-integrable systems, 
such as trapped 87 Rb atoms, which do not noticeably equilibrate ever after thousands 
of collisions [2]. Although the system is in a magnetic trap, this may be a manifestation 
of the quantum integrability of the Lieb-Liniger Hamiltonian for the atoms placed in 
the flat region of the trap. It was firstly conjectured by Rigol et al. |5J that, in dealing 
with integrable systems, one should employ not the usual density matrix (function of 
the energy only) but, instead, the density matrix of a generalized Gibbs ensemble which 
involves all the integrals of motion l\ of the systems 

f> g en~e-^ aih . (1) 

This conjecture has been the subject of several studies, mostly done by using specific 
models [6HTU] , but also studied for initial thermal distributions [11] . An important step 
forward was taken in the work of Barthel and Schollwock [12] in which they generalized 
a previous result by Calabrese and Cardy [13], proving rigourosly that for Gaussian 
initial states and quadratic (fermionic or bosonic) systems the conjecture does hold: 
a finite subsection of an infinite system indeed relaxes to a steady state described 
by a generalized Gibbs ensemble where the extra integrals of motions are simply the 
occupation numbers of each eigenmode. Moreover, they stated (postponing the proof to 
a future paper) that a similar result holds also for Bethe-ansatz solvable models, even if 
in this case the extra integrals of motions appearing in the density matrix do not have 
such a simple physical interpretation. 

In this paper, we tackle the issue of thermalization in integrable systems from a 
different point of view: we focus our attention on continuous systems, i.e. integrable 
field theories in 1+1 dimensions. A reason behind this choice is that the notion of 
quantum integrability can be put on a firm ground only using, as criterium, the elasticity 
of the scattering amplitudes of the excitations of the system [H] and, indeed, this is 
a distinguished property of integrable field theories [15J. Moreover, the rich analytic 
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structure of integrable field theories could help us to obtain results that are harder to 
obtain with other methods. At equilibrium, quantum field theory is a very useful tool 
in order to describe the scaling limit of real world condensed matter systems. If such 
a limit exists also out of equilibrium is still an open question: in fact, after a quantum 
quench a large number of high energy state is populated (HUE], while a field theory 
usually describes only the low energy excitations of the corresponding lattice model. 
However, studying quantum quenches in field theory is still an interesting topic, since 
we expect that the qualitative picture emerging from this analysis should at least give 
precious hints on the behavior of real world systems. The study of quantum quenches 
in integrable systems has also been addressed in [T81IT9] . 

The structure of this article is as follows. In section [2] we briefly remind the main 
properties of integrable field theories and we properly define the kind of quenches we are 
interested in. The main result of this paper is stated in section |3j where we analyze the 
long time limit of the one point function of a local operator. Using a form factor 
expansion, we show that its asymptotic value could be obtained from a physically 
transparent generalized Gibbs ensemble, where the conserved charges are simply the 
occupation number of each eigenmode. Since the argument in section [3] may appear a 
little abstract, in section H] we compute the exact time dependent one point function 
of the energy density operator of the Ising model, showing that its asymptotic value 
agrees with the general formula. Finally, we present our conclusions in section |5j In 



Appendix A| we gather all the technical details of our derivation. 



2. Integrable field theories and quantum quenches 

2.1. Basic facts about integrable quantum field theories 

In this paper we focus our attention on integrable quantum field theories in 1+1 
dimensions (for a review see e.g. [15], [20], [21]). These theories are characterized by an 
infinite set of local conserved currents that greatly constrain the dynamics. Indeed, only 
elastic scattering events can occur in these theories: no particle production is allowed 
and the sets of the final momenta coincide with the initial one. Moreover, due to Yang- 
Baxter equation, the amplitude of the n — > n particle scattering event can be written 
in terms of the 2 — » 2 scattering amplitudes. Using these remarkable properties of the 
scattering processes and some additional assumptions on the analytic structure of the 
S matrix, one is able to determine all the scattering amplitudes and, from their pole 
structure, the mass spectrum of the theory. For sake of clarity in the following we focus 
our attention on the simplest integrable field theories, i.e. those with only one excitation 
of mass m (e.g. Ising model or Sinh-Gordon Model, where the latter is directly relevant 
for one- dimensional bose gases [22l[23]), even if our results could be easily extended to 
any integrable field theory with richer spectrum of excitations. 

A convenient basis can be constructed by the action of a creation operator Z\9) on 
a vacuum state |0), where 9 is the rapidity of the particle (with energy and momentum 
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given by E = mcosh(#), p = msinh(/9) respectively). Since the theory is not free, the 
operators Z(9) and Z\9) do not simply commute or anticommute: they satisfies instead 
the Zamolodchikov-Fadeev algebra 

z{e x ) z{9 2 ) = s(e 1 - e 2 ) z(e 2 ) z(e 1 ) , 

z\e 1 ) z\e 2 ) = s(e l - e 2 ) z\e 2 ) z\e x ) , (2) 
z(e x ) z\e 2 ) = s(6 2 - e 1 ) z\q 2 ) z(9 l ) + 2 u 5(e 1 - e 2 ) , 

where S(9) is the two particles S matrix. Consistency of this algebra is ensured by the 
relationship S(9) = S(—9) and by the unitarity condition S(9)S(—9) = 1. Except for 
the free bosonic theory, where S(0) = 1, the S"-matrix of all other interacting theories 
satisfies 5(0) = —1. Multi-particle states are given by acting with an ordered string of 
the creation operators Z\9i) on a vacuum state |0), i.e. 

\9 1 ,...,9 n ) = tf{9 1 )...tf(9 n )\0), 

(9 n , . . . ,9x\ = <0| Z(9 n ) . . . Z(9x) = (\9 1 ,...9 n ,})^ (3) 

where 

9 1 >9 2 >...>9 n 

These orderings select a set of linearly independent vectors that form a basis in the 
Hilbert space. A convenient resolution of the identity is given by 

'•Ei/s-^ w w 

n=0 J 

where the integration is extended to all values of the rapidities. 

The structure of an integrable theory is so powerful that permits to compute also 
the matrix elements of any local operator 0{x) on the basis ([3]). Indeed, by solving a 
set of monodromy and recursive equations (based on the ^-matrix), one can determine 
the form factors [21] 

F o (9 1 ,...,9 n ) = (0\O(0)\0u---,0 n ). (5) 

Consider, for instance, the Ising field theory at T > T c . This model can be described 
by a free fermionic theory with S'-matrix S = —1 and only one kind of particle of mass 
m oc (T — T c ). There are three physically interesting operators in such a theory: the 
energy operator e (that is proportional to the trace of the energy momentum tensor), 
the spin operator a (the order parameter) and the disorder operator /j, related to a by 
the Kramers- Wannier duality. With a proper choice of the overall normalization, their 
form factors are given by the formulas pi] 

9 | -i27Tm sinh(^f^) ifn = 2 

I otherwise ' 



, {f) p . } °i^' 2 lT< m tanh(^) if „ is odd 

' Mh otherwise ' (?) 

• , H n , j °* n/2 nr< w tanh(^) if n is even 

W>-A)-< otherwise , (8) 
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where a = 25 e~l A 3 mh and A = 1.282427 ... is the Glasher constant) 

It is clear that, from formula (JSJ), we can easily restore the x and t dependance of 
the matrix element (0\O(x, t)\9±, . . . , 9 n ), since the energy and the momentum operators 
are diagonal in the basis ([3]). However, it is not so obvious how to obtain from dSJ) the 
generic matrix element (9 n , . . . , 9i\O(0)\9[, . . . , 9' n ): this task can be accomplished by 
exploiting the crossing symmetry. If A and B are two sets of rapidity, we have 

(A\O(0)\B) = S AAl S BBl (At\O(0)\Bi)(A 2 \B 2 ), (9) 

A=AiUA 2 ;B=BiUB 2 

where the sum is over all the possible ways of splitting the sets A/B in two subsets A\j 
Bi and A 2 /B 2 while S^Ai an d Sbbx are the products of S(9) we need to rearrange the 
rapidities in the proper order, namely 

(A\ = S AAl (A2A 1 \, (10) 
\B) = SbbABA) . (11) 

The symbol A± in denotes that each rapidity 9\. . .9 r in A\ is shifted by an 
infinitesimal amount 6j so that (A^\0\Bi) is simply related to the form factors (jSJ) 

(^+|O(0)|B!) = (0\O(0)\BiA+ -iit). (12) 

Here comes however the tricky point. If the e, are finite, the form factors are (for 
real rapidities) regular functions, while if we take the 6j — > limit (as we need to 
do at the end of the calculations) the form factors usually diverge, because we are in 
the kinematical situation in which some of the rapidities of the bra and the ket states 
coincide: the simplest case of this circumstance is provided by the 2-particle matrix 
element (9\fi\9') where \x is the disorder operator of the Ising model (jHl), which indeed 
diverges when 9 = 9'. This discussion shows that a prescription is needed for handling 
these kinematical divergencies. The one proposed in [251126] consists of taking only the 
regular part of f|T2|) and discarding all the terms proportional to an inverse power of q 

(9 n ,...,9 1 \O(0)\9[,...,9' m ) conn = (13) 



Finite Parts 



lim (O|0(O)|0i, ...,9' m ,9 n -i7r + e n ,...,9 1 -i7T + e 1 



It should be stressed, however, that this prescription alone is not enough to properly take 
care of all the divergencies and, usually, it must be supplemented with extra corrective 
factors coming from the Bethe-ansatz technique [25j. Without entering in too many 
details, to show this aspect, consider for instance as a density matrix of the syste nfl 

= exp{-fd9\(9)h(9)) 
Z 

where h(9) = Z^(9)Z(9) and X(9) is an appropriate function of 9. Notice that, if 
A(#) = 7p m cosh(#) the density matrix (114j) describes the familiar canonical ensemble, 

f Here we are using the conformal normalization of the operators, such as in the limit x — > 

(<f>(x,0)4>(0,0)) ~ r^ps where A is the conformal dimension of the field. 

| [n(6) , n(0')] = 0, so we do not have ordering problems in the exponential. 
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while, in the more general case, it can be associated to the generalized Gibbs ensemble 
dl]): indeed, in an integrable field theory the (infinite) set of local continuity equations 
determines the conserved charges Q s labeled by an integer index s [15] , and each charge 
Q s can be written as 



Q s 



d9e sd h{9), 



(15) 



where q s is a constant. Obviously, it would be nice to write down the average of a local 
operator O w.r.t. the density matrix ( TH|) in term of these form factors. But, if we do 
it applying the prescription ([131 alone we end up with the result 

(O)p = Tr(pO(0)) = ! WRONG! (16) 



+oo 

£ 

71=0 



1 f d9 1 ... d9 r , 



nl 



(2*0 



n 

i=l 



-A(0i) 



1 - S{0) e-W) 



! WRONG ! 

,0110(0)10!,.. 



. 9 n )c 



which is wrong since it does not agree with the thermodynamic Bethe ansatz. It was 
firstly conjectured by LeClair and Mussardo |25j that the correct expression is instead 



(O)* = Tr(pO(0)) 



(17) 



£ 

n=0 



1 f d9 l ... d9 r , 



nl 



n 

i=l 



-A(0i) 



1 - S(0) e-^) 



9 1 \O(0)\9 1 ,...,9 n ) 



i w n/ conn i 



(19) 



where the A are dressed according to the integral equation 
\(9) = X(9) - J^<P(0- 9') log[l + e"^')], 
ip being the derivative of the phase shift 

Originally, ( JT7|) was only an educated guess but its correctness was confirmed by 
subsequent checks [27-29J. Even if presently there is still no complete proof of the 
exactness of ffl71) . a step forward its rigorous derivation has been taken in the recent 
papers by Pozsgay and Takacs [30H32] . who showed that ffl7|) might be obtained in a 
rigorous way by putting the system in a finite volume, thus discretizing the rapidities 
to regularize the kinematical singularities, and taking the infinite volume limit only at 
the end of the calculations. 



2.2. Integrable boundary condition and quantum quenches 

Our aim is to study the time evolution of the expectation value of a local operator O(x) 
on a pure state \B) that is not an eigenstate of the integrable Hamiltonian H, i. e. 

«?(»■«)).= <MEp g> . m 

As shown by Calabrese and Cardy [13, 33], performing a Wick rotation, this dynamical 
problem can be mapped into a statistical problem defined in a strip geometry, where 
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the initial state \B) plays the role of a boundary condition. In an integrable field theory, 
the most natural boundary states are the ones that do not spoil the integral of motions. 
These integrable boundary states were originally studied by Ghoshal and Zamolodchikov 
in [33] and their properties can be summarized as follows: these states have a peculiar 
form, given by a coherent superposition of pairs of equal and opposite rapidity (Cooper 
pairsjj 



\B) = exp 



1 r+oo in 

- / ^-K"»{9)Zl{-e) Z\{9)) 

2 ./ ^ 2 IT 



|0>, (21) 



where the indices a and b run over all the particles. The amplitude K ab {6) satisfies a set 
of equations that depends on the S matrix (such as boundary Yang Baxter, boundary 
unitarity and crossing equations): the different solutions of these equations provide the 
set of integrable boundary conditions of the theory in question. As an explicit example, 
let's consider once again the Ising model, in which there are three integrable boundary 
conditions associated to the amplitudes 

Kfi X {6) = ztanh(|) Fixed Boundary Condition 

Kfree{9) = — i coth(|) Free Boundary Condition (22) 

Kh{0) = i tanh(|) j^oshle) Magnetic Boundary Condition, 

where k = 1 — ijf- and h is a real number (the boundary magnetic field). Clearly, the 
magnetic boundary condition interpolates between the free and the fixed ones. So, we 
see that only particular choices of the amplitude K give us an integrable boundary state. 

Our interest in the following is to study a particular class of quenches, i.e. those 
where the initial state is expressed by a coherent superposition of particle pairs, as 
the one given in eqn.f l2~Tj) . Among these quenches there are, of course, the integrable 
boundary states but the whole class of these quenches is actually larger than the genuine 
integrable ones: indeed, for theories with only one particle, it suffices that the amplitude 
K{9) satisfies the "boundary cross- unitarity equation"[jj] 

K{9) = S(29)K(-9), . (23) 

This equation simply comes from the request of invariance of the expression (1211) under 
a change of variable 9 — > —9. 

Notice that, at least for a free fermionic theory as the Ising model, a boundary 
state of the form ( |2T|) has a very simple meaning within the theory of quench processes. 
Suppose in fact that %o is a free fermionic theory, with its ground state \B) identified 
by the condition Z (9)\B) = 0, where Z (9) is the annihilation operator related to 
Hq. If, at t — 0, the Hamiltonian is suddenly changed from Ho to a new almost free 
Hamiltonian TL (for instance, making a quench of the mass from m to m), the new 
sets of creation/annihilation operators Z\Z will be related to the initial ones by a 
Bogoliubov transformation of the form 

Z (9)~ [Z{8)-f{-9)Z\-9)\, (24) 

§ Here, for simplicity, we ignore the presence of additional zero-rapidity terms in the expression of \B). 
|| For theories with many particles, K a f,{9) must also solve an appropriate boundary Yang-Baxter 
equation to ensure commutativity of the various terms in the exponential expression (|21[) . 
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Hence, the algebraic equation that identifies the boundary state in terms of the new 
creation/annihilation operators is given by 

Z (6) \B) = =^ [Z{9) - f{-9)Z\-9)] \B) = (25) 

whose solution is 



\B) = exp 



d9 
2^ 



f(9)Z\-9) tf(0)) 



|0>. 



(26) 



Summarizing, in the following we consider those quantum quenches where the initial 
state is given by a superposition of Cooper pairs according tolSTl). where the amplitude 
K is a (regular) function that satisfy (123 p . However the boundary state as written in 
(|2ip is typically not normalizable. Indeed, if the amplitude K does not goes to zero 
for large rapidities (as in eq.( l22]) ). it means that we are exciting modes with arbitrary 
high energy, hence the state has an unphysical divergent energy density. So, as done 
for the critical systems [TBI 133] , it is therefore convenient to introduce an extrapolation 
time To that make the norm of the state (12~T1) finite: this parameter plays the role of an 
ultraviolet cut-off, in any case present in physical systems. So, we are interested in 

(B\e int O(x)e- iHt \B) 



where 



(0{x,t)) 1 



\B) = e- Hro \B) 



(B\B) 



exp 



d9 
2^ 



G(9)Z\-6) Z\9)) 



|0>, 



(27) 



(28) 



and G(9) = e 2mr ° cosh( . e ) x \9) . Notice that the extrapolation time does not spoil the 
basic equation of the boundary state, namely also the new amplitude G(9) satisfies (1231) 

G{9) = S(29) G{-9) and G{9) = S(-29) G(-9) . (29) 



3. Long time behavior 

In this section, we are concerned with the computation of 

0= lim (0{x,t))s , (30) 

t— !-+oo 

and in the following we will show that O can be expressed as an average over a density 
matrix like ( TH|) . At first sight, it seems very unlikely that we can do so: our boundary 
state (128]) has a very peculiar structure, since it is the superposition of pairs of opposite 
rapidity, while in an average over an ensemble (1171) there is no sign of such a structure. 
How comes that the system retains no memory of this pair structure in the long time 
limit? 

First of all, a trivial remark: since \B) is translational invariant, (0(x,t))^ do not 
depend on x and therefore, from now on, we will set x = 0. In principle, it is quite clear 
what we have to do in order to compute (1301) : first we expand the exponential in ( 1281) 
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and, taking into account that the Hamiltonian is diagonal in the particle basis, we thus 
arrive to the double sum 



(0{x,t)) 



+00 

J^y — 

Rl R\ n\ 1 



d9 x ... dd n d6[ . . . de[ 



(4 



7T 



(4tt) 



exp [2it(E n (9)-E l (9'))] 



11 G ^ 



i=l 



U=i 



h, —0i\O\ - 9 1 ,9' 1 , ... — 9{, 6[) 



(31) 



where we used the short-hand 



E n {9) 



m 



n 

E 

i=i 



cosh(6'j) . 



(32) 



However it is difficult to compute the long time limit directly from fl31~]) . The reason is 
that the matrix element in ( l31~i) are not regular functions (rather they have delta-like 
contributions) and, in such a case, we cannot apply a stationary phase argument. In 
order to isolate the singular parts, we have to employ eqn.(j2]). Consider, for instance, 
the term with n — I — 1 in the numerator of fl3~D). i.e. 



df) dff 

exp [2it{E x {9) - E 1 (9'))]G(9 1 )G(9' 1 )(9 1 , -9 t \0\ - 9[,9' 1 ) . (33) 

4 7T 4 7T 



Applying the crossing relation ()9]), we can recast this term as 
dQx d9[ 

4 7T 47T 

d9 x 



e^p[2it(E l (9) - E l (9'))]G(9 l )G(9[)(9t,-9t\0\ - 9',,9'x) + (34) 



2?r 



\G(9 1 )\ 2 (9+\O\9 1 ) + (0\O\0) 



An An 



A) 



where we make use of the symmetry properties f)29p of G. Actually, the inner product 
(9i,— 9± I — 9'i,9[) above is divergent in the infinite volume limit and it should be 
regularized by putting the system in a box of length L, obtaining 

m L 



(5(9-9')Y 



2tt 



cosh(9)5(9-9'). 



(35) 



This, however, is not an important contribution since, as shown in Appendix A, all 
these inner products cancel out with the denominator of (I3TI) . What is really crucial 
is that, apart from these infinite volume divergencies that we can easily regularize, the 
integrands in (1341) are all well-behaved functions: hence we can now easily take the 
infinite time limit t — > +00, so that fl34l) simply becomes 



2n 



|G(^)| 2 (^|O|^) + (0|O|0) 



^f^G(9 1 )G(9' 1 ){9 1 ,-9 1 \-9[,9' 1 ) (36) 

4 7T 4 7T 



because the first term in (1341) vanishes for the fast oscillation of its integrand. 

In the light of this example, the strategy to compute the expectation values of local 
operators can be stated as follows. 

(i) We first expand the exponential in the numerator of ( I28p . ending up with the double 
sum ASH). 
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(ii) Then, we use in order to isolate the delta-like terms and, after having done that, 
we take the infinite time limit, where all terms that explicitly depend on time go 
to zero, due to the fast oscillation of the integrand. This is a simple consequence 
of the stationary phase argument, that can also be seen in the following way. If 
the infinite time limit exists (and the stationary phase argument assures us that it 
does exist), then it must coincide with the temporal average 

= lim (0(x,t)) 6 = lim i ^ dt'{0(x,t')) s . (37) 



We know that the time-dependent part of the numerator of (1311) consists of a sum 
of terms as 

J dd x ... dd n d0' x . . . d9\ exp [2 % t{E n (9) - E0))] F{9 1 ...9 n ,9[... 9[) , (38) 

where F is the regular function obtained by applying So, the only contributions 
to the infinite time limit comes from the region E n {9) = Ei(9'), whose Lebesgue 
measure is zero, so the integral goes to zero since F has no delta-like term. 

With these steps in mind, it is a nice combinatorial exercise to show that O can be 
finally expressed as 



O 



+oo 1 „ n 



n=0 



d9< 



i (27T) 



\G{6i 



Ll-S(0)|G(^)| 2 J 



9 n +ie n ...0 1 +ie 1 |O(O)|0 1 ...0 n }.(39) 



In order to keep our exposition clear, we will present the details of the combinatorics 
behind (l39|) in the Appendix A However, ( )39|) is still a meaningless expression, since 
we have to regularize it in a proper way. One way to do it is by analogy with LeClair 
and Mussardo formula, discussed in section 12.11 When we perform an average over a 
density matrix (I14p . we end up with the following expression 



n=0 



I'll 



TT — 



-\(6) 



1 



9 n +ie n ...9 1 +ie 1 \O(0)\9 1 ...9 n ).{40) 



LeClair and Mussardo suggested that the proper way to regularize the — > limit of 
this expression is to take the the connected part of the form factors ( T131) and to dress 
X(9) according to the integral equation fflBl . This regularization scheme holds for every 
function A(6>). The situation is the same in equation ( 1391) . with |G(#)| 2 that plays the 
role of e~ x( > 9 \ So, the natural way of regularize ( )39|) lead us to 

12 



o 



r 1 



n=0 



1 r d9j 



\G{9 % 



...9 l \O(0)\9 1 ...9 n ) ( 



(41) 



S(0)|G(^)| 2 . 

where |G(6 I )| 2 is dressed in the same way as the term e~ x ^ entering the thermodynamic 
Bethe ansatz 

dff 



\G(9)\=\G(9)\ A ^ 



2tt 



ip{9-9') \og[l + \G{9)[ 



(42) 



The above dressing formula is based on the LeClair and Mussardo conjecture (it has 
actually the same mathematical structure) and on the possibility to exchange the — > 
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and t — > +00 limits. Assuming that such a regularization scheme is indeed correct, it 
turns out that that the long time limit of ff27|) could be described by a generalized 
Gibbs ensemble f fT4"l) . where the constant of motions are simply given by the occupation 
number h{6) and the function X(9) is fixed by the conditions 

(*!W=T rfc »(«)), (43) 
(B\B) 

thus proving Rigol et al.'s conjecture for integrable field theory. If we look at our starting 
point f )27|) . this result is quite unexpected: when we expanded the exponential in ff27|) 
we had a double sum, and it was only thanks to the infinite time limit that we could 
rewrite it as a single summation. Moreover, while the boundary state (I28I) is formed 
by pairs of particles with opposite rapidity, this feature is completely lost in the final 
expression ( 14~TT) . 

4. The simplest example 

It is instructive to see how the general ideas of the previous section apply in the simplest 
case provided by the one point function of the e operator of the Ising model ([6]). Indeed, 
for this operator we can calculate exactly its one point function for any time with 
elementary techniques. From the form factors (j6j), it follows that the operator e is a 
quadratic form in the creation - annhilation operators 

e (0) = /^{2 7 rmcosh(^) *(A)Z(A) + 

[Z{^)Z{h)+Z\h)Z\^)\ ). (44) 



% tc m sinn 



+ — ■ 2 

Since the theory is free, we can easily calculate the expectation value of binomials of 
the creation-annihilation operators on the boundary states (introducing, for instance, a 
generating functional) and we have that 

(B(i)|Zt(ft)Z(ft)|B(t)> |G(ft)| 2 

m = (A " ■ 2> i+igw (45) 

<SjS) = 2 * 4 < ft + ft) l + |G(A)|» ' (46) 



Hence 



rf6> |G(fl)| ; 
2^1 + |G(^)| 2 

c/6» sinh(^) Jm exp (-2 z 



<«(,,«)>,= + 2, m y- r L^L_ + (47) 
— 2 7rm 



2vr l + l^)! 2 

The results (j4Tj) has all the features we expect to hold in the general case. First of 
all, the time dependent part goes to zero as a consequence of the fast oscillation of the 
integrand. The long time asymptotic value obviously agrees with our general result: 
in this case, the structure of the operator is so simple that the entire sum reduces 
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to a single term. More important, since we are able to calculate exactly the time 
dependence, we can also estimate the approach to the t — > +00 limit value using 
a stationary phase approximation. It turns out that ( 1471) approaches its asymptotic 
value as an inverse power law, in contrast with the exponential decay of the massless 
case (33] [13]. Heuristically, the different behavior of the two cases (massive and 
massless) can be understood in the following way. We know that the relaxation to 
the stationary value is due to the speed of the quasiparticle excitations: in a massless 
theory, all the quasiparticles move at the fastest possible velocity, so we can expect that 
the decay to the limit value is the fastest possible (exponential decay), while, in a massive 
theory, the time decay has to be slower for the spread of the velocity distribution of the 
quasi-particles. It may be interesting to compare this continuum formula f l47j) with the 
one point function of the transverse magnetization of the quantum Ising model [35J: 
a quench of the transverse magnetization in the lattice corresponds to a mass quench 
in the continuum. In both cases there are oscillations modulated by a power law with 
exponent — |. However, it must be noted that the long time leading behavior of the 
transverse magnetization on the lattice is not only dictated by the small momenta, and 
therefore it is not a surprise that the continuum theory is unable to capture all the 
features of the lattice model. 

5. Conclusions 

In this paper we studied the quantum quenches in the context of integrable field theories 
in 1+1 dimension. We considered a specific class of quenches, those whose initial state 
is a coherent superposition of Cooper pairs (|28|) . a structure that is in agreement with 
the integrability of the theory in the bulk. For such quenches, we analyzed the one point 
function of a local operator and argued that, in the long time limit, it can be neatly 
expressed in terms of the formulas (jjl] ) and (HSl) : the final result employs a series over 
the connected form factors ( fl~3l) . which are finite expressions. 

As discussed extensively in the text, our final result must be presently considered 
only as a solid conjecture since it is based on a regular izat ion scheme of the kinematical 
singularities of the original form factors. This regularization is closely related to the one 
used in LeClair and Mussardo formula for the one point function at finite temperature 
[25] - formula that has successfully passed many checks although it is not yet rigorously 
proved. In this respect, we have performed an additional but independent test of our 
result for a particular one point function of the Ising model. If the full correctness of 
eqn. (fill) will be confirmed by a future rigorous analysis, it follows that the asymptotic 
value of the one point function of a local operator could be computed as an average over 
a density matrix ( fl4l) with a temperature different for each eigenmode: in short, our 
result is a strong evidence that Rigol et a/.'s conjecture does hold for integrable field 
theories. 

In this paper we limited ourselves to integrable theory with only one massive particle 
but our results could be extended to theories with more than one particle without any 
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physical difficulty, although with a combinatorics much more involved. Finally, our 
results seems to suggest that in a generic massive integrable theory the decay towards the 
asymptotic value is dictated by a power law (in contrast with the conformal exponential 
decay) whose leading exponent can be determined by a stationary phase approximation. 
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Appendix A. A diagrammatic proof of eq. 



In this appendix we would like to show that, for t — > +00, 
(B\0(x,t)\B) 



(B\B) 

£^1 fd9 1 



n=0 



d9 r . 



(2tt)™ 



n 



i=i 



\G(0i 



(A.l) 



(e:...et\o(o)\e 1 ...e n ) 



1-S(0)\G(9^\ 

where 9^ = 6{ + zej. In this appendix, we will call the form factors like 
(9£ ■ ■ ■ 9f\O(0)\9i . . . 9 m ) regular form factors, in order to distinguish them from the 
complete form factors (9 n . . . 9x\O(0)\9i . . . 9 m ), that have also delta-like contributions. 
For €j finite these regular form factors are continuous functions but, despite their name, 
they can have a singular — > limit, hence the need of the regularization procedure 
previously discussed. Let's firstly briefly anticipate the main steps of the proof: after 
we expand the numerator of the l.h.s. of (1A.1I) (as done in (|3~T|) ). we obtain 



J'=l 



?n 3 — On, ■ ■ ■ Q\> —9\\0\ — 1 ,9 V . . . — 9[, 9[) . (A. 2) 



Then, applying repeatedly the crossing relation (Q to the matrix elements in (1A.2I) . we 
arrive to an expression in which we can take the infinite time limit. In this limit, (1A.2I) 
reduces to the r.h.s. of flA.lj) times the denominator of the l.h.s 



n=0 



where 



Z = (B\B) = (B\e iUt e- iHt \B) = J^ Z n , 

r 

d9 x ... d9 n d9[ . . . d9' n 



Z " n ~ 12 
nr 



\\G{9 l )G{0[) 



(4vr) n (47r) n 

• exp(2 % t (E n (9) - E n (9'))(9 n , -9 n , . . .6 1 , -9 X \ - 9[, 9[, . . . - 9' n , 9' n ) 
thus completing the proof. 



(A.3) 



(A.4) 
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Figure Al. In this figure we present the building blocks of our diagrammatic 
representation. On the left, we can see a bra and a ket with n = 3 Cooper pairs. 
In the middle, we see one of the diagrams that represent the inner product of the bra 
and the ket. The figure on the right shows a term proportional to the regular form 
factor (02 , &i\O\0i, 0' 2 )- The circle in the middle stands for the operator O, that has 
2r = 4 legs, half connected to the bra and half to the ket. 



Introductory remarks 

Here we would like to highlight some basic facts that we will find useful for our proof. 
First of all, we notice that, when we recast ( 1A.2I) in terms of the regular form factors, 



the only terms that survive in the infinite time limit are the time-independent ones. As 
a consequence, we can discard all the terms in the double sum in ( 1A.2I) where n ^ I. 
Moreover, from (Q it is clear that the term (0|(9|0) in the r.h.s. of flA.lj) is correct, so 
in the following we will not consider this contribution. 

In order to follow more easily our ideas, it may be useful to develop a diagrammatic 
representation. In fig. IA1I it is shown a bra state with n pairs (a ket state can 
be introduced in a similar way). An inner product between a bra and a ket both 
made of n pairs can be represented as the sum of all the possible ways to link 
together the particles of the ket with the particles of the bra, each link meaning a 
contraction between the corresponding particles. Of course, one should be careful about 
the permutation of particles and the corresponding S matrix, but we will take care 
later of these details. We are interested in the matrix elements of the operator O 
between states made of Cooper pairs. In particular, our aim is to reduce the full 
matrix element (9 n , —9 n , . . . 9±, — 9\\0\ — 9[,9' 1 , . . . — 9' n , 9' n ) to the regular form factors 
(9f r , ... ,9^1019'^, ... ,9 j r ) . These regular terms can be diagrammatically represented as 
in fig. IA11 where the operator O has 2r legs: r of them are connected to the particles in 
the bra, while the other r are connected to the particles in the ket. So, our combinatorial 
problem reduces to a problem in which we have to connect the legs of the operator to 
the bra/ket and the remaining particles together, according to ([9]). In order to get the 
right combinatorial coefficients, we have to remember what follows. 

• As we stated before, only the matrix elements with the same numbers of particles 
in the ket and in the bra give a non vanishing contribution in the long times limit. 
Therefor, the number of legs connected to the bra is always equal to the number 
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of legs connected to the ket and, from now on, we will only specify the number of 
particles connected to the bra. 

• In principle, when we consider the matrix element between two states with n pairs, 
i. e. (9 n , —9 n , . . .9\, —0i\O\ — 9[, 9' x , . . . — 9' n , 9' n ), we could expect to end up with a 
sum of regular terms with r legs linked to the bra, with r < 2n. However, it turns 
out that r < n: if we connect r particles to the operator, we are left with 2n — r 
delta functions, and to suppress all the time dependencies, we have to eliminate n 
integration variables, hence r < n. 

• Finally, if we link a particle to the operator, its pair partner cannot be connected 
to O, otherwise their time dependence survives. 



The disconnected terms 

In this subsection we will show how some contributions to flA.lj) (the disconnected terms) 
cancel out with the denominator Z ( 1A.3I) . In order to understand this point, we analyze 
the matrix elements {62, —92, #1, — 9i\0\ — 9[, 9[, —9' 2 , 9' 2 ) . As it is shown in fig. \K2\ we 
have essentially two types of diagram. Let's focus our attention on the second one: it is 
clear that the contribution in the dotted box factorize from the integral with the form 
factors. Hence, if we have a set of particles that is completely disconnected from the 
operator O (this means that no one of these particles or their partners is connected to 
or to a particle whose partner is connected to O), its contribution factorizes. 

What we want to show now is that these disconnected pieces cancel out with the 
denominator Z. Let us consider the term in flA.2j) with n pairs in the bra (let's denote 



it as O n ), and let's focus our attention on the contribution 0[ n> k\, such that k pairs in 
the bra (as well as k pairs in the ket) are disconnected and so only n — k pairs in the bra 
are connected to O. Of course, a similar term can also be obtained from O n ^k, when 
no particle is disconnected. With our notation, this term can be written as O[ n _k,o] ■ If 
we are able to show that, for any n and k < n, 

0[ ny k] = C[n-fc,0] %k , (A.5) 
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then it follows that, when we sum over all n and k < n, (1A.2j) becomes 

'+00 \ 
\m=0 / 

hence the disconnected pieces cancel out with Z. 

The proof of (1A.5[) is actually quite simple. It is clear that the inner product of 



the k disconnected pairs gives an integral proportional to Zj~, as well as 0[ n ,fe] is clearly 
proportional to O[ n _k,o]- So, in order to complete the proof, we have only to check the 
proportionality constant. 0[ n} k] has an overall coefficient from the expansion of the 

exponentials, while we can choose the k disconnected pairs in ^ fc ,^ fc ^ equivalent ways, 
since the creation operator of a pair commutes ([Z^ (—9) Z^ (9) , Z^ (—9') Z^ (9')] = 0). 
Nicely, ^ n \y2 is the correct coefficient for O[ n _k,o] and is the right one for Z kl thus 
concluding the proof of IA.5I 

The last step 

Here we conclude our proof, showing that 0[ m> o] has actually the right structure to give 
f lA.ip . Let's call 0{ m ^ the contribution where the operator has r legs connected to 
the bra, the bra has m pairs and no particle is disconnected from O. We have already 
pointed out that for t — > +00 

in 

C[m,0] = ^C{m,r}, (A-7) 
r=0 

since the sum is restricted to r < m. So, in order to get ( 1A.1I) . we need only to show 
that 

o {m , r} = -<! dd ) n -f r {et...et\o(o)\e 1 ...9 r }- (A.8) 



r! J (2vr) 

• ' ^(or^ac^)! 2 ) 41 . . . s(o) v - i (|G(6> r )i 2 ) ir ] 



ti,...,i r 



where the summation J is over all the positive integers ij such that Y2j ij = m - 

In order to prove (1A.8|) . we need to be a little careful with the ordering of the particles 
and the labeling of the rapidities. However, if we exchange two particles, the contribution 
is same, since (as we already know) the pairs do commute while the exchange of two 
particles forming a pair is equivalent to change of the integrable variable 8 — > —9. This 
is a consequence of the symmetry of G fl29|) . 

Before concluding our proof of ( 1A.1I) . we need to understand how to label the 
rapidity. We start with 2m integration variables and the delta-functions reduce them to 
r. So, we use the convention to label the rapidities as in (1A.8j) : we call 9\ the rapidity 



of the particle of the bra closest to O and we take advantage of delta functions in such 
a way that the rapidity of the particle in the ket nearest to O is also #1, and so on. 

We can now finally show that we obtain exactly the structure ( 1A.8I) . First of all, 
from our previous reasoning about the long time limit, it is clear that, performing all the 
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possible contractions, for the term in the second row of (I A . 8 [) we arrive to an expression 



as 

' [c^xicwrr ...crWOGwit] > (a.9) 

ii,...,i r 

where ci(ii) . . . c r (i r ) are unknown constant. What we have to prove is that 

(i) The overall coefficient agrees with ( 1A.8I) 

(ii) <•,(/,) sw ! . 

The first point comes out from the combinatorial coefficient that takes in account all 
the equivalent way to link the particles and the operator. In O[ m) o] , we have an overall 
coefficient that is ^2^: the factorial comes from the exponentials while we get the | 
from the integration measure that is ^- and not We can choose the r particles in 
the bra (and the r in the ket) connected to O in 

J 2 (A ' 10) 

ways. We remind that if a particle is connected to O, its pair companion cannot be 
directly connect to the operator, otherwise the contribution is time dependent hence it 
goes to zero for long times. We have to determine in how many ways we can connect the 
particle in the bra with rapidity 61 to the others in order to have a term like (|G(#i)| 2 ) n . 
The answer is 

[2(m — r) 2(m — r — l)...2(m-r-i! + l)] 2 r, (A.ll) 

where the r comes out from the r equivalent ways to choose a particle in the ket 
connected to O. If we repeat the same argument for all the particles, we end up with 
an overall coefficient that is 

1 1 9 2r 1 

^_ m ,2 2 2(m-r) r , = 1 ? 



2m2(m — 1) . . . 2(m — r — 1) 



m \2 r \2 r \ 

in agreement with IA.8I 

Finally, we show that Cj(ij) = S'(O)^ -1 . Since we know that every permutation of 
particles gives the same contribution, it is sufficient to show it only for one of the many 
equivalent ways to link particles. In particular, we will consider the following way to 
separate the rapidities in two sets 

{9 m , —9 m , ■ ■ ■ ,0\, —6i\0\ — 9[, 9[, . . . — 9' m , 9 m , ) = Saa 1 Sbb 1 (9 r , . . . 9 1 \0\9' 1 , . . . 9' r ) ■ 
(9 m , —9 m , ■ ■ ■ , 9 r +i, —9 r +i, —9 r , . . . — 9\ | — 9[, ... — 9 rl —9 r+1 , 9 r+1 , . . . , — 9 m , 9 m ) . (A. 13) 

An useful trick is to remember that the contractions are such to have, at the end, 
9i = 9\ for i = 1, . . .r. When we impose this condition, we see that the S matrices in 
(IA.I3j) reduces to the identity. Now, we want to contract —9\ with —9' r+l , obtaining a 



5(9i — 9' r+1 ) and the desidered IG^i)! 2 . Then, we link —9 r+ i to 9' r+1 . This contraction 
gives us a ^(^i + 9 r+ \). Finally, we commute 9 r+ i (that now is equal to —Oi) with 
— 9 r . . . — 9 2 . In this way, we end up in a situation similar to the initial one. We have 
two pair less, a |G(#i)| 2 overall and a huge product of S matrices that comes from all 
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the exchanges done. However, if we remind that at the end of our calculation we have 
9i — 6[ for i — 1, . . .r, it is easy to see that this product of S matrices reduces to S(0) 
. We can repeat this procedure until the overall coefficient is (S'(0)|G(^i)| 2 ) il . Then 
we contract —Q\ with —9[, obtaining another (G^x)! 2 and we go on doing the same 
manipulations on —62. It is clear that at the end we obtain exactly fA.81 
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